/* Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010 Mario Orsi
   This file is part of Brahms.
   Brahms is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
   Brahms is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
   You should have received a copy of the GNU General Public License along with Brahms.  If not, see <http://www.gnu.org/licenses/>. */

// utils.c -- various useful functions

#include "dataStructs.h"

/* From Rapaport [pp.491-492]: */
#define IADD   453806245
#define IMUL   314159269
#define MASK   2147483647
#define SCALE  0.4656612873e-9

extern int moreCycles;

int randSeedP=17;

real **CreateRealMatrix( int n1, int n2 ) 
{
  int i, j; /* counters */
  real **a; /* a pointer to a vector of vectors ( = a matrix). */

  /* Ask for memory space for the array of pointers to arrays. 
     Notice the sizeof(real *) is not the size of a real, but the
     size of a pointer to a real. */
  a = (real **)malloc(n1 * sizeof(real *));

  /* Now allocate all memory required for the whole matrix. Let 
     the returned pointer be the pointer to our first array. */
  a[0] = (real *) malloc (n1 * n2 * sizeof(real));

  /* Now for the funny part. Calculate the rest of the needed pointers
     ourselves. The pointer to the second array (your second row or
     column) must be placed n2 elements ahead of the first array. This
     is untrivial pointer arithmetic. */
  for ( i = 1; i < n1; i++ )  {
    a[i] = a[i-1] + n2;
  }

  /* Set all elements in the matrix to zero. */
  for ( i = 0; i < n1; i++)  {
    for ( j = 0; j < n2; j++) {
      a[i][j] = 0.0; /* remember: 0 is an integer, 0.0 is a float */
    }
  }

  /* Return with the address of the array of arrays: */
  return a;
}

void InitRand( int randSeedI ) /* Rapaport [p.492] */
{
  struct timeval tv;

  if ( randSeedI != 0 ) randSeedP = randSeedI;
  else {
    gettimeofday( &tv, 0 );
    randSeedP = tv.tv_usec;
  }
}

void InvertMatrix( RMat *I, const RMat *A ) // does what it says but hey: it assumes to operate with a 3 x 3 matrix!
{
  real d = A->u[0]*(A->u[4]*A->u[8]-A->u[7]*A->u[5])-A->u[3]*(A->u[1]*A->u[8]-A->u[7]*A->u[2])+A->u[6]*(A->u[1]*A->u[5]-A->u[4]*A->u[2]);
  I->u[ 0 ] = ( A->u[ 4 ] * A->u[ 8 ] - A->u[ 7 ] * A->u[ 5 ] ) / d;
  I->u[ 3 ] = ( A->u[ 6 ] * A->u[ 5 ] - A->u[ 3 ] * A->u[ 8 ] ) / d;
  I->u[ 6 ] = ( A->u[ 3 ] * A->u[ 7 ] - A->u[ 6 ] * A->u[ 4 ] ) / d;
  I->u[ 1 ] = ( A->u[ 7 ] * A->u[ 2 ] - A->u[ 1 ] * A->u[ 8 ] ) / d;
  I->u[ 4 ] = ( A->u[ 0 ] * A->u[ 8 ] - A->u[ 6 ] * A->u[ 2 ] ) / d;
  I->u[ 7 ] = ( A->u[ 6 ] * A->u[ 1 ] - A->u[ 0 ] * A->u[ 7 ] ) / d;
  I->u[ 2 ] = ( A->u[ 1 ] * A->u[ 5 ] - A->u[ 4 ] * A->u[ 2 ] ) / d;
  I->u[ 5 ] = ( A->u[ 3 ] * A->u[ 2 ] - A->u[ 0 ] * A->u[ 5 ] ) / d;
  I->u[ 8 ] = ( A->u[ 0 ] * A->u[ 4 ] - A->u[ 3 ] * A->u[ 1 ] ) / d;
}

real RandR() // yields a random value in the range (0,1). [Rapaport, p.492].
{
  randSeedP = ( randSeedP * IMUL + IADD ) & MASK;
  return ( randSeedP * SCALE );
}

static void ProcInterrupt() // allows orderly termination via, eg, $kill -s SIGUSR1 {PID} [Rapaport, p.494] 
{ 
  printf("Orderly termination via user signal.\n");
  moreCycles = 0; 
}

void SetupInterrupt() // Rapaport [p.494] 
{ 
  signal( SIGUSR1, ProcInterrupt ); 
}

void VRand( VecR *p ) /* Rapaport [p.492] */
{
  real s, x, y;

  s = 2.;
  while ( s > 1. ) {
    x = 2. * RandR() - 1.;
    y = 2. * RandR() - 1.;
    s = Sqr( x ) + Sqr( y );
  }
  p->z = 1. - 2. * s;
  s = 2. * sqrt( 1. - s );
  p->x = s * x;
  p->y = s * y;
}
